{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "93b34000",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import random\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "195e187f",
   "metadata": {},
   "outputs": [],
   "source": [
    "D = {}\n",
    "path = ''\n",
    "with open(path+'OpenAlexID_Year_DindexAllYears.txt','r') as f:\n",
    "    for line in f:\n",
    "        line = line.strip('\\n').split('\\t')\n",
    "        D[line[0]] = float(line[-1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "c439cafd",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "49417708"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(D)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "a4308105",
   "metadata": {},
   "outputs": [],
   "source": [
    "D_sample = random.sample(list(D.values()),1000000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "d49b7b36",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAucAAADoCAYAAABMzko0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAApqklEQVR4nO3df1yV9f3/8efxB2I/9KuijFSEyq9pmCYnWuhKzTBRaCZkP0RZS3Nsu1F+a+vjdwtzlvtsn3TOho5WYmtL+qGtY6zCSUrqapD4zdpWVgaCxg9TY4koXN8/zBMHkB/n13Wdw+N+u3HjnOt6n+t6+vbiOi/evK/r2AzDMAQAAADAdD3MDgAAAADgLIpzAAAAwCIozgEAAACLoDgHAAAALILiHAAAALAIinMAAADAInqZHcBKwsLCFBUVZXYMAAAABKiDBw+qpqbG7ddTnDcTFRWl4uJis2MAAAAgQNntdo9eT3EuyeFwyOFw6Pjx42ZHAQAAQDfGnHNJSUlJysnJUf/+/c2OAgAAgG6M4hwAAACwCIpzAAAAwCIozgEAAACL4IJQAECHpv/iVZfnr/98pklJACC4UZyLu7UAAADAGpjWIu7WAgAAAGugOAcAAAAsguIcAAAAsIigLM7j4+M1btw4xcTEaPny5WbHAQAAADolKC8Ife2119SvXz81NjZq0qRJSk5O1vjx482OBQAAALQrKEfO+/XrJ0lqaGhQQ0ODyWkAAACAzrFccb5z504lJydr6NChstlsys3NbdUmOztb0dHRCg0NVWxsrIqKilq1ufbaazVkyBBNmzaNUXMAAAAEBMsV53V1dYqJidGaNWvUt2/fVuvz8vKUmZmppUuXau/evYqPj9eMGTNUVlbm0u7tt99WRUWFSktLtX//fn/FBwAAANxmueI8MTFRjz32mFJSUtSjR+t4q1atUnp6uhYuXKjRo0dr7dq1ioiI0Lp161q17devn6ZOnarXXnvNH9EBAAAAj1iuOG9PQ0ODSkpKlJCQ4LI8ISFBu3fvliQdO3ZMNTU1kqT6+nq98cYbuuKKK867zZycHNntdtntdlVXV/suPAAAANCBgLpbS01NjRobGxUeHu6yPDw8XNu2bZMkHT16VKmpqTp9+rSampp02223adasWefd5qJFixQRESGHw6HCwkKf5gcAAADaE1DF+Tk2m83luWEYzmWXXnqpSkpKurS9pKQkJSUlyW63ey0jAAAA0FUBVZyHhYWpZ8+eOnLkiMvyqqqqVqPpXeFwOORwOHT8+HFPIwIAAABuC6g55yEhIYqNjVVBQYHL8oKCAsXHx7u93aSkJOXk5Kh///6eRgQAAADcZrmR87q6Oh04cECS1NTUpLKyMpWWlmrgwIGKjIzUkiVLlJaWpri4OE2cOFHr169XZWWlFi9e7PY+GTkHgNbu/t2bZkcAgG7HcsV5cXGxpkyZ4nyelZWlrKwsLViwQLm5uZo7d65qa2u1YsUKHT58WDExMcrPz9eIESPc3idzzgEAAGAFlivOJ0+eLMMw2m2TkZGhjIwMr+2TkXMAAABYQUDNOfcV5pwDAADACiw3cm4GRs4BQFr8+51mRwCAbo+RczFyDgAAAGtg5BwAurEf/+EtsyMAAJqhOBfTWgAAAGANFOfiVooAuo/7NuwyOwIAoB3MOQcAAAAsgpFzAAhyS3J3mx0BANBJFOdizjmA4PLgM3t8vo8HNn6zj/9ZcJ3P9wcA3QXFuZhzDiDw/fTZv5sdAQDgBRTnAACPPPTs2y7PfznvWpOSAEDgozgHgAC09E9vd9wIABBwKM4BIED83z+/Y3aETvn5c9/k/MUdcSYmAYDAQ3EOABb1803/MDsCAMDPKM7F3VoAmGdZXrHzsWEYJibxjaxm/z5JemQuF94DQHtsRjC+G7jJbreruLi444YA4KZHnnc9xzQ/A7c8HTe1eG1Tk3Heto3NnjdvJ0mNzZ43tXxd4zfPTze67rHhzDfPK47+x2Xd2MiBzse9erp+nl3vnjbn4x49XNc1W6WHb6NQBxB8PK0nGTkHAC/7xQslLs8ZAQEAdBbFOQDAFCtedP0l5mcpsSYlAQDrCLrivLy8XGlpaaqqqlLv3r2VlZWlW2+91exYAILAys17z7uOGYKea96//3Xr1SYmAQDzBF1x3qtXL/3mN7/R+PHjVVVVpdjYWN1888264IILzI4GIMD8csv5i3H4Vsu+f2g2xTqA7iHoivOIiAhFRERIkoYMGaIBAwaopqZGkZGRJicDYEW/ernU5Tkj4NbU/P/pJ98db1oOAPC1Hh038a+dO3cqOTlZQ4cOlc1mU25ubqs22dnZio6OVmhoqGJjY1VUVNTmtoqLi3X69GkNHz7cx6kBBJJfvVzq/ELg+Z+/7HP5AoBgYrmR87q6OsXExGj+/PmaP39+q/V5eXnKzMxUdna2Jk2apOzsbM2YMUMffPCBy+h4bW2t5s+fr6eeeko2m63VdgAEn8dfabtQYyw8uLX8f/8/yeNMSgIAnrP0fc4vuugiPfHEE0pPT3cuu/baa3XVVVfpySefdC4bOXKkUlJStHLlSknSqVOndNNNN2nhwoVKS0vr9P64zzkQeFZv/X9nHxit7+F9TsulLe8D7tK22TZatmrvbNn8dY0ttm+cp13LbXa3+5z36OE6cNKzWduWYyo9mi3o0WJle23VYt39s64SAPhSt7rPeUNDg0pKSvTAAw+4LE9ISNDu3bslnX3DSk9P19SpUztVmOfk5CgnJ0eSVF1d7f3QALxuzavvSWJ+OLru3LEjSZkzx5qYBADaFlDFeU1NjRobGxUeHu6yPDw8XNu2bZMk7dq1S3l5ebrqqqv08ssvS5L++Mc/auzYtk/CixYtUkREhBwOhwoLC32aH4Bnfpv/XseNgE5qXqhLFOsArCGgivNzWs4hNwzDuWzSpElqamr5x+D2JSUlKSkpSXY7HyUNWM0Tf91vdgR0E+0daz+aEePHJAC6s4AqzsPCwtSzZ08dOXLEZXlVVVWr0fSucDgccjgcOn78uKcRAXjButffl2FIBpdywiLWvf7+edf9YPqVfkwCINhZ7laK7QkJCVFsbKwKCgpclhcUFCg+Pt7t7SYlJSknJ0f9+/f3NCIAN61/4wNJ7RdBgBWtf+MD5xcAeMpyI+d1dXU6cOCAJKmpqUllZWUqLS3VwIEDFRkZqSVLligtLU1xcXGaOHGi1q9fr8rKSi1evNjtfTJyDpjr918XNb+nuEGAa3UMN5uFee9NY/wbBkBAstytFN98801NmTKl1fIFCxY4P5AoOztbv/rVr3T48GHFxMRo9erVuv766z3eN7dSBPznD9v+KensbQabms5OYLHpm1v9nZvW0vwM1Xzd2e/NJr5wK8U227XM091vpdjyta77d29d8+ugWrVq/u9tmbvZ43umjT7v9gEEFk/rScsV52Y4N3JeWFiojz76yOw4QNB66m//cinAJYpzivNm67pxcd4SxToQuLrVfc59hbu1AL719PZ/tVuIAHDV/Gfme1OvMDULAP/qUnFeUFCgm266yVdZTMOcc8A3cgv/TVEOeOjcz9GCKaPMjgLAD7o0raVHjx669NJLde+99yo9PV2DBw/2ZTa/Y8454Lk/7vhQ0tnpE4Zx9k/3jYbx9ZQVMa2FaS3N1jGt5Xx69LA519tsXz+2nd2GYUjzJ//vdl4NwEye1pNdupXi9u3bdc011+jnP/+5hg8frjvvvFM7duxwe+cAgsezOz/Ss18X5gB869kdH+qPOz7kZw4IQl2a1jJ58mRNnjxZtbW1evrpp/Xkk09q06ZNGjVqlBYvXqz58+drwIABvsoKwIL+XMRF1ICZ/rTzm59Bm0268zsjTUwDwFNufQjRoEGD9OCDD+rDDz9UQUGBwsLCtGTJEg0dOlTp6el67733vJ3TpxwOhxYtWsScc6ALNr11QM+9dcDsGABa4BdmILB59Amh+fn5+u1vf6u///3vGjJkiObPn68dO3ZowoQJWrdunbcy+hyfEAp07PndHzsf5+2iKAesbNPXvzjn7TqgvF0fd9AagJV0uTg/cuSIHn30UUVHR2vWrFk6duyYnn32WZWXl2v9+vU6cOCA7r33Xi1fvtwXeQGY4IWvC/MXdn/sUqQDsK7mP6sv8HMLBIwuzTmfM2eOtm7dqtDQUM2bN08ZGRm68sorXdr07NlTd955p7Kzs70aFID/vfT3T8yOAMBLXtjzsfOOMXO+fale+vsnmvPtS01OBaClLhXnH330kX7zm98oLS1NF1100XnbjR07VoWFhR6H8xfucw642kxRDgS1zW9/anYEAOfRpeJ869atioiIUO/evVutO3PmjCorKxUZGamLL75YN9xwg9dC+hqfEAp84+V3eNMGupO/fP0zf0tctMlJAEhdnHMeHR2tvXv3trlu3759io7mBxsIVK/846D+8o+DZscA4Ed/afbL+Cv/OKhXOAcAputScd7eh4mePn261SfBAQgMjuLPzI4AwCIcxQfNjgB0ax1Oazl27JiOHj3qfF5RUaFPPnGdj3ry5Elt3LhR3/rWt7yfEIBPvFryWauPpweAc7aWfOa8gDRxQqTy3y1T4oRIk1MBwa/D4nzNmjV65JFHZLPZZLPZlJKS0mY7wzD0yCOPeD2gP3BBKLqT/HfL9PX7LQC0Kf/dsjaX/3VvmWZcTYEO+FKHxfl3v/tdRUVFyTAM3X333frZz36myy67zKVNnz59NGbMGF111VU+C+pLXBCKYPZ6abkkafr44c7HANAVr+1tXay/Xlqu6eOHm5AGCG4dFufjxo3TuHHjJEk2m00zZ85UWFiYz4MB8Nwb+74pxgv2HTIxCYBgUbDvkJq+vgatYN8h3TRumMmJgODSpSs4FyxYQGEOAACctv2/Qy7fAXimw5HzqVOnKjs7W1dccYWmTp3ablubzaa//e1vXgvnruTkZBUVFenGG2/Uiy++aHYcwO+2v1dhdgQAAOCGDkfOm98+sampSYZhnPerqanJp2E76/7779czzzxjdgzAb97cXylJKtxfocL9FOYA/OvcgMCb+yud56M33680MxIQsDocOS8sLHQ+fvPNN32ZxWumTJkSMFkBb+GNEIBVND8f7Xi/UjdceYnzO4D2We5Tg3bu3Knk5GQNHTpUNptNubm5rdpkZ2crOjpaoaGhio2NVVFRkf+DAiYr+udhl+8AYDVFH3B+ArqqS8X5X/7yF23YsMH5/LPPPtN1112niy++WCkpKaqrq/M4UF1dnWJiYrRmzRr17du31fq8vDxlZmZq6dKl2rt3r+Lj4zVjxgyVlbV9T1YgmOz61xGX529RmAMIIJyzgI51qThfsWKFqqurnc+XLFmiQ4cOadGiRdq5c6eWLVvmcaDExEQ99thjSklJUY8ereOtWrVK6enpWrhwoUaPHq21a9cqIiJC69at83jfQKDY/e8jHTcCAIvY8+/PnY+bn7/OPW6+Hq6aDD7LubvpUnH+8ccfOz9o6OTJk8rPz9eqVav0+OOP67HHHtOWLVt8EvKchoYGlZSUKCEhwWV5QkKCdu/e7dY2c3JyZLfbZbfbXX7xAKzg7x9+3uZjAAh05wpyzm2Aqy4V5/X19c6pJrt379aZM2echfKoUaNUWenbC9JqamrU2Nio8PBwl+Xh4eE6cuSb38SnTZum1NRU5efna9iwYdqzZ895t7lo0SJlZWVpwoQJCgkJ8Vl2wFve/og3MgCBrbPnMc536I66VJxHRUXprbfeknR2/nlsbKz69+8vSaqqqnI+9jWbzeby3DAMl2Xbtm1TdXW1vvrqKx06dEjXXXddu9tLSkpSTk6O3/IDLb1XdtTl+T8OVJmUBADMUfxx5/563dl2QKDq8FaKzd1777164IEHtGXLFpWWlrrM896zZ4/GjBnj9YDNhYWFqWfPni6j5NLZXwxajqZ3hcPhkMPh0PHjxz2NCLit9GCNxkeFqeQT3ngAdE8ln1TLJpsmXBqmkjaK8LaWAcGmSyPnmZmZys3N1XXXXaenn35aCxcudK778ssvlZ6e7u18LkJCQhQbG6uCggKX5QUFBYqPj3d7u4ycwyzvlx8977rST2ucj9/9pOa87QAgmO39lPMfupcujZxL0l133aW77rqr1fLf//73XglUV1enAwcOSDr7iaRlZWUqLS3VwIEDFRkZqSVLligtLU1xcXGaOHGi1q9fr8rKSi1evNjtfTJyDqvYd7C21bLSg7wxAUBz+w7WalzUoFaPgWDQ5eL8nKqqKtXX17daHhkZ6VGg4uJiTZkyxfk8KytLWVlZWrBggXJzczV37lzV1tZqxYoVOnz4sGJiYpSfn68RI0a4vc+kpCQlJSXJbrd7lB3orH8e+kKjhw0wOwYABKz3PqvV2BFni/L3yo5KhuF8/vGR47rsW23/Nbyi9j+SpKGDLvRPUKCLulScnzhxQpmZmcrLy9OpU6fabNPY2OhRoMmTJ8vo4J6eGRkZysjI8Gg/zTFyDl87cPi4Lo9g2hQAAGhfl4rzH/7wh3rppZf0/e9/X2PHjlWfPn18lcuvGDmHJypq/+MyAlNWU6fIsIs6fN0/D33h8h0AcFZ71+MAwa5Lxfnrr7+uX//61/rhD3/oqzxAUDhY9aWihlxsdgwACBr/qmAgA91Dl4pzwzA0atQoX2UxDdNaAAAAYAVdupXi7bffLofD4asspuFWivCnjw7zSyAAdAXnTXQnXRo5T0hI0H333acvv/xSiYmJGjhwYKs2U6dO9Vo4wNvqG84oNMTtmxR57OMjvMEAgLe0dU79+MgJE5IA3tOlKuWWW26RJH366afKzc11LrfZbDIMQzabzeO7tZiBaS0AAACwgi4V54WFhb7KYSru1gJf+Kz6S40YfLEOVn0pSc7vAADf4DyLYNCl4vyGG27wVQ4gKH1WzRsFAJihs7e1BaymSxeEnlNTU6OtW7dq48aNOnr07L1I6+vr1dTU5NVwAAAA7jr3aaAtHwNW1qXi3DAMPfjggxo2bJiSk5N199136+DBg5LOzkd/9NFHfZHR5xwOhxYtWsSccwAA0ErlUQp7+E+XivOVK1fqiSee0MMPP6y3335bhmE41yUlJWnr1q1eD+gP3EoR3nL4i6/MjgAAAAJYl+ac/+EPf9DDDz+s//qv/2p1V5bLL79cH3/8sVfDAYHiyDGKcgCwos+PfaXw/3WBJKnq+MkO21cdP6kh/ft6PYfZt/JF4OjSyHlFRYW+/e1vt7kuJCRE//kPf/YBAAAA3NWl4nzo0KHav39/m+v27dun6Ohor4QCAkHNiXrVfllvdgwAgJs4h8OKulScp6amavny5dq1a5dzmc1m04cffqjHH39ct99+u9cD+gMXhAIAELy+qDvlle1QzMMfulScL1u2TFdccYWuv/56jRw5UtLZgn3s2LEaOXKkHnroIZ+E9DUuCEVHOjNPEQAAwFNdKs779u2rN998Uxs3blR8fLymTZuma665Rjk5OSooKFBISIivcgKWcOJkw9nvXzWYnAQA4Ilz53FPzue8F8AXunTZcH19vYqLi9WnTx9997vfVUREhGJjYxUaGuqrfIDl1NWfNjsCAMAHvjz5zfn9xMkG9evLoCP8r1Mj56dOnVJmZqYGDhyoG264Qbfffrvmzp2r66+/XoMGDdIDDzyghgbr/PaYn5+vUaNGaeTIkcrOzjY7DgAAANApHY6cG4ahWbNmafv27brllluUmJioyMhIGYah8vJybd26VatXr9YHH3yg/Px8f2Ru15kzZ5SZmant27dr0KBBstvtmj17tiIiIsyOhgDT8p609Q1nTEwDAAC6gw6L8xdffFGFhYV68cUXNXv27Fbr77nnHm3evFm33XabNm/erFtvvdUnQTvrnXfe0ZgxYzR8+HBJ0uzZs7V161YtXLjQ1FwAACCwnDrdqD69e5odA91Mh9NannvuOd12221tFubn3HrrrUpNTdWf/vQnjwPt3LlTycnJGjp0qGw2m3Jzc1u1yc7OVnR0tEJDQxUbG6uioiLnusrKSmdhLknDhg1TRUWFx7nQfTWcaey4EQCgWzjd2NRhG9434IkOi/O9e/dq5syZHW5o1qxZevfddz0OVFdXp5iYGK1Zs0Z9+7b++Ny8vDxlZmZq6dKl2rt3r+Lj4zVjxgyVlZVJOjsNpyWbzeZxLnRPZzpxEgYABL7m53vO/TBTh8V5dXW1IiMjO9xQZGSkqqqqPA6UmJioxx57TCkpKerRo3W8VatWKT09XQsXLtTo0aO1du1aRUREaN26dZLOfoppeXm5s/2hQ4d0ySWXeJwLAADAHU1tDBwC59Nhcf7VV1+pT58+HW4oJCRE9fW+/eSshoYGlZSUKCEhwWV5QkKCdu/eLUmKi4vT+++/r/Lycp08eVJbtmxpd+Q/JydHdrtddrtd1dXVPs0PAAAAtKdT9zmvqKjQJ5980m6bQ4cOeSVQe2pqatTY2Kjw8HCX5eHh4dq2bZskqVevXlq9erVuvPFGNTU1KTMzs92R80WLFikiIkIOh0OFhYU+zQ8AAAC0p1PFeUpKSodtDMPw29zulvtpue+kpCQlJSV1envn2tvtdq9lBAAAcNfpxib17tmlD3JHkOiwON+wYYM/cnRKWFiYevbsqSNHjrgsr6qqajWa3hUOh0MOh0PHjx/3NCL85Exjk3px0gIA4LyaDEM9uClGwOmwOF+wYIE/cnRKSEiIYmNjVVBQoNTUVOfygoICzZkzx+3tMnIOAAAAK+jUtBZ/qqur04EDByRJTU1NKisrU2lpqQYOHKjIyEgtWbJEaWlpiouL08SJE7V+/XpVVlZq8eLFbu+TkXMAAABYgeWK8+LiYk2ZMsX5PCsrS1lZWVqwYIFyc3M1d+5c1dbWasWKFTp8+LBiYmKUn5+vESNGuL1PRs4BAABgBZYrzidPntzmBwk1l5GRoYyMDK/tk5Hz4MLHLQMAgEDFFXU6O3Kek5Oj/v37mx0FAADAp/hQJGujOAcAAAAswnLTWszAtBYAABAMGpua1LMHY6+BjP89Ma3FbB1dYwAAANzH+2xgoTgHAAAALIJpLWJaCwAAAKyBkXMxrQUAAADWQHEOAAAAWATFOQAAADzCRafew5xzMeccAAAA1sDIuZhzDgAAAGugOAcAAAAsguIcAAAAsAiKcwAAAMAiuCBUXBAKAAAAa2DkXFwQCgAAAGugOAcAAOjmGpu6fp9y7m3uGxTnAAAAgEUEbXGenJysAQMGKCUlxewoaKaJ37IBAADOK2iL8/vvv1/PPPOM2TEAAACATgva4nzKlCm6+OKLzY4BAAAAdJrfi/OdO3cqOTlZQ4cOlc1mU25ubqs22dnZio6OVmhoqGJjY1VUVOTvmAAAAIDf+b04r6urU0xMjNasWaO+ffu2Wp+Xl6fMzEwtXbpUe/fuVXx8vGbMmKGysjJnm5iYmDa/ysvL/flPAQAAALzK7x9ClJiYqMTERElSenp6q/WrVq1Senq6Fi5cKElau3atXnvtNa1bt04rV66UJO3fv99veQEAAAB/sdSc84aGBpWUlCghIcFleUJCgnbv3u2Tfebk5Mhut8tut6u6uton+wAAAAA6w1LFeU1NjRobGxUeHu6yPDw8XEeOHOnStqZNm6bU1FTl5+dr2LBh2rNnT5vtFi1apKysLE2YMEEhISFuZwcAAAA85fdpLZ1hs9lcnhuG0WpZR7Zt29bptklJSUpKSpLdbu/SPgAAAABvslRxHhYWpp49e7YaJa+qqmo1mu5NDodDDodDx48f99k+AAAAgI5YalpLSEiIYmNjVVBQ4LK8oKBA8fHxPttvUlKScnJy1L9/f5/tAwAAAOiI30fO6+rqdODAAUlSU1OTysrKVFpaqoEDByoyMlJLlixRWlqa4uLiNHHiRK1fv16VlZVavHixzzIxcg4AAAAr8HtxXlxcrClTpjifZ2VlKSsrSwsWLFBubq7mzp2r2tparVixQocPH1ZMTIzy8/M1YsQIn2VizjkAAACswO/F+eTJk2UYRrttMjIylJGR4adEjJwDAADAGiw159wszDkHAACAFVCcAwAAABZhqVspmoVpLQAAALACRs7FtBYAAABYA8U5AAAAYBFMaxHTWgAAAGANjJyLaS0AAACwBopzAAAAwCIozgEAAACLYM65mHMOAAAAa2DkXMw5BwAAgDVQnAMAAAAWQXEOAAAAWATFOQAAAGARXBAqLggFAACANTByLi4IBQAAgDVQnAMAAAAWQXEOAAAAWATFOQAAAGARNsMwDLNDWEVYWJiioqIkSdXV1Ro8eLC5gQIUfecZ+s999J3v0ce+Rf96hv5zH33nmeb9d/DgQdXU1Li9LYrz87Db7SouLjY7RkCi7zxD/7mPvvM9+ti36F/P0H/uo+88483+Y1oLAAAAYBEU5wAAAIBFUJyfx6JFi8yOELDoO8/Qf+6j73yPPvYt+tcz9J/76DvPeLP/mHMOAAAAWAQj5wAAAIBFUJwDAAAAFkFxDgAAAFgExbmX5Ofna9SoURo5cqSys7PNjhPwkpOTNWDAAKWkpJgdJaCUl5dr8uTJGjNmjMaNG6fNmzebHSngxMfHa9y4cYqJidHy5cvNjmMZ2dnZio6OVmhoqGJjY1VUVGR2pG6PY9V9nCs9w3t057lTH3JBqBecOXNGo0eP1vbt2zVo0CDZ7Xb97W9/U0REhNnRAlZhYaHq6uq0ceNGvfjii2bHCRiHDx/W559/rvHjx6uqqkqxsbH697//rQsuuMDsaAHjxIkT6tevnxobGzVp0iStW7dO48ePNzuWqfLy8jRv3jxlZ2dr0qRJys7O1oYNG/TBBx8oMjLS7HjdFseq+zhXeob36M5xtz5k5NwL3nnnHY0ZM0bDhw/XBRdcoNmzZ2vr1q1mxwpoU6ZM0cUXX2x2jIATERHhfHMeMmSIBgwY4NFHCHdH/fr1kyQ1NDSooaHB5DTWsGrVKqWnp2vhwoUaPXq01q5dq4iICK1bt87saN0ax6r7OFd6hvfoznG3Puz2xfnOnTuVnJysoUOHymazKTc3t8127f1Jt7KyUsOHD3c+HzZsmCoqKnwd3RTe6K/uyt99V1xcrNOnT7scm4HMn/137bXXasiQIZo2bVq3H4lsaGhQSUmJEhISXJYnJCRo9+7dJqWyts4cq976OQ/GY9Wf/ScF17nS330X7DztT3frw25fnNfV1SkmJkZr1qxR375922yTl5enzMxMLV26VHv37lV8fLxmzJihsrIySVJbM4NsNptPc5vFG/0lSTExMW1+lZeX++uf4nf+7Lva2lrNnz9fTz31VNAci/7sv7ffflsVFRUqLS3V/v37ff5vs7Kamho1NjYqPDzcZXl4eLiOHDliUipr6+hY9eY5MhiPVX/2X7CdK/3Zd92Bp/3pdn1owOnCCy80NmzY0Gp5XFyccc8997gsu/zyy42HHnrIMAzD2LVrl5GcnOxct3TpUiMnJ8enWa3A3f7qrMLCQmPOnDmeRLQsX/ZdfX298Z3vfMd45plnPI1pWb4+9s755S9/afz6179267XBoqKiwpBk7Ny502X5smXLjFGjRpmUKnC0dax6+zg1jOA9Vn3Zf8F+rvT1sRfM79Ftcac/3a0Pu/3IeUc68yfduLg4vf/++yovL9fJkye1ZcsWzZw504y4puNP4O7zRt8ZhqH09HRNnTpVaWlpvohpWd7ov2PHjjnnndbX1+uNN97QFVdc4fWsgSQsLEw9e/ZsNUpeVVXVajQdHfPWObK7Hqve6r/ueK7k/dm7fFkfUpx3oDN/0u3Vq5dWr16tG2+8UWPHjtUPfvADXXLJJWbENZ23/gQ+bdo0paamKj8/X8OGDdOePXu8HdVyvNF3u3btUl5enl5++WWNHz9e48eP13vvveeLuJbjjf47evSopk+frquuukp2u1033HCDZs2a5Yu4ASMkJESxsbEqKChwWV5QUKD4+HiTUgUub50ju+ux6q3+647nSm9OUeuO79Et+bI+7OWTxCb72c9+pkcffbTdNoWFhZo8eXKnt9lyjpBhGC7LkpKSlJSU1KWcVmFGf3Vk27ZtnW5rJqv13aRJk9TU1NTpfZnNav136aWXqqSkpNP76i6WLFmitLQ0xcXFaeLEiVq/fr0qKyu1ePFis6MFLE/Pkd39WPW0/wLtXOlNnvadFDjv0f7gi/owKIvz++67T/PmzWu3TWfvzdsd/qRLf7mPvvMM/RcY5s6dq9raWq1YsUKHDx9WTEyM8vPzNWLECLOjBRyOU8/Qf+6j77zLl/0ZlMV5WFiYwsLCvLKt5n/STU1NdS4vKCjQnDlzvLIPs9Ff7qPvPEP/BY6MjAxlZGSYHSPgcZx6hv5zH33nXb7sz6Aszruirq5OBw4ckCQ1NTWprKxMpaWlGjhwoHPEjj/pfoP+ch995xn6D4Gio2OV47R99J/76DvvMq0/3b2lTLAoLCw0JLX6WrBggUu73/3ud8aIESOMkJAQY8KECcaOHTvMCWwy+st99J1n6D8Eis4cqxyn50f/uY++8y6z+tNmGG3cIR0AAACA33ErRQAAAMAiKM4BAAAAi6A4BwAAACyC4hwAAACwCIpzAAAAwCIozgEAAACLoDgHAAAALILiHAACUG5urmw2m/PrwgsvVFRUlGbPnq3nn39eTU1Nnd7GwYMHvZYrKipK6enpXtseAHQ3vcwOAABw3wsvvKBhw4bp1KlTKisr06uvvqo77rhDOTk5cjgc6tu373lfO3PmTO3Zs0cRERF+TAwAaA/FOQAEsPHjx+vyyy93Pk9LS1NqaqpSU1P1k5/8RGvXrj3vawcPHqzBgwf7IyYAoJOY1gIAQWbOnDm65ZZb9OSTT+qrr746b7u2prVERUVp3rx52rRpk0aPHq0LL7xQdrtdb731VqvXr1mzRlFRUQoNDZXdbldRUVGb+/n000911113afDgwerTp4/Gjx+vLVu2ONfv379fffv21f333+/yuqVLl6pPnz7au3dvF3sAAAIXxTkABKHExESdOnVKxcXFXX5tUVGRHn/8cf3iF79QXl6eGhsbNWvWLB07dszZ5qmnntJ9992nKVOm6OWXX1Z6erruuOMOffHFFy7bKi8v17XXXqt9+/Zp9erVeuWVVzRhwgTNmTNHr7zyiiQpJiZGjz/+uNasWaO//vWvkqTCwkL993//t1auXKmrr77a/Y4AgADDtBYACEKRkZGSpMOHD3f5tSdOnFBpaakGDBggSfrWt76la665Rvn5+brzzjvV1NSkZcuWafr06dqwYYPzdYMHD9btt9/usq1ly5bJMAzt2LFDgwYNkiRNnz5d5eXlevjhh5WcnCxJysjI0BtvvKH09HRt375d8+bN00033dRqNB0Agh0j5wAQhAzDkCTZbDY1NjbqzJkzzq9z687nuuuucxbmkjR27FhJUllZmSTp0KFDOnTokG677TaX182ZM0e9ermO+bz22mtKTExU//79XTJMnz5d+/bt04kTJ5xtn3rqKfXu3VuxsbE6c+aMNm7cKJvN5n4nAEAAojgHgCBUXl4uSYqIiNBll12m3r17O782btzY7msHDhzo8rxPnz6SpPr6eknfjMaHh4e7tOvVq5dzdPycqqoqPfPMMy777927tx588EFJUm1trbPtoEGDNHPmTJ06dUp33HFHq+0DQHfAtBYACEKvvvqqQkNDFRsbK4fDoVOnTjnXRUdHe7Ttc7de/Pzzz12WnzlzxqXYls4W3N/5znf005/+tM1tXXLJJc7H27Zt05NPPim73a7s7GzNmzdPdrvdo6wAEGgozgEgyGzevFmvvPKKMjMzdcEFFzinpXjLsGHDNHz4cD3//PO6++67nctfeuklnTlzxqXtzTffrD179ujKK69s957rNTU1mj9/vhITE7VlyxZNmjRJd955p959911ddNFFXs0PAFZGcQ4AAay0tFQ1NTVqaGhQWVmZtm7dqhdeeEE33XSTVq5c6ZN99ujRQ1lZWbrnnnv0ve99T7fffrsOHDiglStXql+/fi5tly9frri4OF1//fX60Y9+pKioKH3xxRfav3+/PvnkEz399NOSpLvvvluGYWjDhg3q3bu3/vznP+vqq6/Wj3/8Y5eLTgEg2FGcA0AAS01NlSSFhoZqyJAhmjBhgjZt2qSUlBSfXkz5/e9/X3V1dVq1apWee+45xcTEaNOmTZo3b55Lu8jISBUXF2vZsmVaunSpqqurNWjQIMXExGjBggWSpCeeeEJbt27V66+/7vxQpMsuu0zZ2dlKS0vTzTffrLlz5/rs3wIAVmIzOrpsHwAAAIBfcLcWAAAAwCIozgEAAACLoDgHAAAALILiHAAAALAIinMAAADAIijOAQAAAIugOAcAAAAsguIcAAAAsAiKcwAAAMAi/j8i3TTy09GERgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x216 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "nbins=50000\n",
    "fig = plt.figure(figsize=(10, 3),facecolor='white')\n",
    "ax = fig.add_subplot(111)\n",
    "#\n",
    "N, bins, patches = ax.hist(D_sample,density=True,bins=nbins,edgecolor='white', linewidth=0)\n",
    "LN=np.log(max(N))\n",
    "\n",
    "\n",
    "# Find fractions of small natural denominator and make them less distractive\n",
    "delta=bins[1]-bins[0]\n",
    "peaks=set([j/i for i in range(1,100) for j in range(1,i)]+\\\n",
    "[-j/i for i in range(1,100) for j in range(1,i)]) \n",
    "ys=dict((i,(patches[i-1]._height+patches[i+1]._height)*0.5) \\\n",
    "        for i in range(1,len(N)-2)) #smooth color\n",
    "\n",
    "for i,z in zip(range(0,nbins),bins):\n",
    "    if patches[i]._height>0: # only plot patches with values\n",
    "        r=(np.log(patches[i]._height)/LN)**1 # transparency scales to y values\n",
    "        for p in peaks:\n",
    "            if z<p<z+delta and 0<i<len(N):\n",
    "                r=(np.log(ys[i]+1)/LN)**1 # disruption fraction\n",
    "        r = max(0, min(1, r))\n",
    "        patches[i].set_alpha(r)\n",
    "        patches[i].set_facecolor(color='SteelBlue')\n",
    "\n",
    "#\n",
    "plt.yscale('log')\n",
    "plt.xscale('symlog',linthresh=0.001)\n",
    "plt.tight_layout()\n",
    "plt.xlabel('D-index',fontsize=16)\n",
    "plt.ylabel('Density',fontsize=16)\n",
    "plt.xticks(size=14)\n",
    "plt.yticks(size=14)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a272fbc7",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
